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ABSTRACT 



Context. In 2D-simulations of thin gaseous disks with embedded planets or self-gravity the gravitational potential needs to be 

smoothed to avoid singularities in the numerical evaluation of the gravitational potential or force. The softening prescription used 

in 2D needs to be adjusted properly to correctly resemble the realistic case of vertically extended 3D disks. 

Aims. We analyze the embedded planet and the self-gravity case and provide a method to evaluate the required smoothing in 2D 

simulations of thin disks. 

Methods. Starting from the averaged hydrodynamic equations and using a vertically isothermal disk model, we calculate the force to 

be used in 2D simulations. We compared our results to the often used Plummer form of the potential, which runs as oc l/(r^ -I- e^)''^. 

For that purpose we computed the required smoothing length e as a function of distance r to the planet or to a disk element within a 

self-gravitating disk. 

Results. We find that for longer distances e is determined solely by the vertical disk thickness H. For the planet case we find that 

outside r x H & value of 6 = Q.IH describes the averaged force very well, while in the self-gravitating disk the value needs to be 

higher, e = \.2H. For shorter distances the smoothing needs to be reduced significantly. Comparing torque densities of 3D and 2D 

simulations we show that the modification to the vertical density stratification as induced by an embedded planet needs to be taken 

into account to obtain agreeing results. 

Conclusions. It is very important to use the correct value of e in 2D simulations to obtain a realistic outcome. In disk fragmentation 

simulations the choice of e can determine whether a disk will fragment or not. Because a wrong smoothing length can change even the 

direction of migration, it is very important to include the effect of the planet on the local scale height in 2D planet-disk simulations. 

We provide an approximate and fast method for this purpose that agrees very well with full 3D simulations. 
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1. Introduction 

Numerical simulations of accretion disks are often performed 
in the two-dimensional (2D) thin-disk approximation, because 
a full 3D treatment with high resolution is still a computation- 
ally demanding endeavor requiring a lot of patience. In numer- 
ical disk models with embedded planets and/or self-gravity the 
gravitational potential (force) needs to be smoothed because it 
diverges for very short mutual distances. For that purpose the po- 
tential is softened to avoid the singularities that arise, for exam- 
ple, from pointlike objects, such as planets embedded in disks. 
In full 3D simulations this smoothing may be required solely for 
stability purposes and can be chosen to be as small as the given 
numerical resolution allows. In contrast, 2D disk simulations are 
typically based on a vertical averaging procedure that leads to a 
physically required smoothing. Ideally, this smoothing should be 
calculated in a way that the 2D simulations mimic the 3D case as 
closely as possible. The most often used potential smoothing has 
a Plummer form with T oc -l/(r^ + e^)'/^, where r is the distance 
to the gravitating object and e is a suitably chosen smoothing or 
softening length. If a vertically stratified disk has a thickness H, 
we would expect that somehow e should depend suitably on H. 
In a sense, the potential is 'diluted' in this case due to the disk's 
finite thickness. 

For the planet-disk problem, iMivoshi et al.l (Il999h have 
shown with local shearing sheet simulations owing due to this 
dilution effect, the total torque exerted on a planet in a three 



dimensional disk is only about 43% of the torque obtained in 
the thin, unsmoothed 2D case. These authors also showed that 
the strength of the one-sided torque depends on the valu e of the 
smoot hing, where larger e lead to smaller torques. Later. Massej 
(l2002h has studied the smoothing problem in greater detail. He 
has shown that good agreement of the total 2D and 3D Lindblad 
torques can be obtained for smoothing lengths of e = 0.75H, 
whe re H is t he vertical scale height of the disk, see Eq. ( fTTt be- 
low. iMassetl found that for the Lindblad torques this optimum 
e/H value is independent of the planet mass and the thickness 
of the disk. On the contrary, for the corotation torques that are 
generated by material moving on horseshoe orbits, he found that 
the required smoothing depe nds on th e ratio Rh/H, where R^ is 
the Hill radius of the planet. iMassetl concluded that there is no 
'magic' value of e that generates overall agreement of 2D with 
3D results. 

Based on these studies, in 2D planet-disk simulations the 
parameter e is typically chosen such that the total torque act- 
ing on the planet, which determines the important migration 
speed, is approximately equal to that obtain e d thro ugh 3D (lin- 
ear) analysis, for example by iTanaka et al.l (l2002h . This argu- 
ment has led to the choice of e ^ 0-3 - Q.6H, a value very of- 
ten u s ed in these simulations (Masset 20021 : Ide Val-Borro et al.l 
l2006t iPaardekooper & Papaloizoij|l2009h . However, recently a 
very small smoothing has be en advocated for 2D planet-disk 
simulation (iDonget al.ll2011h . A smoothing based on a verti- 
cal integration using Gaussian density profiles has been used 
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by iLi et al.l (l2005l l2009h . but they provided tio details oti the 
methodology and accuracy. 

For self-gravitating disks, the conditions for fragmenta- 
tion have recently attracted much attention in the context of 
planet formation via gravitational insta bility. Typical l y, stud- 
ies are performed in full 3D (see e.g. iMeru & Batel (1201 Ibb 
and references therein). However, to save computational effort, 
2D simulations present an interesting alternative in this context 
dPaardekooper et al1l201 ll) . Here, the accurate treatment of self- 
gravity is very important. 

The incorporation of self-gravity in 2D t hin-disk calculations 
can be achieved by fast Fourier transforms. iBaruteau & Massea 
(l2008h presented a method where the force is calculated di- 
rectly using a smoothing length that scales linearly with radius, 
and they used t he same smoothing, e = 0.3//, for the planet 
and self-gravity. iLi et al.l (12005) used this method to calculate 
the potential in the disk's midplane. The simultaneous treat- 
ment of an embedded planet and disk-self-gravity can be im- 
portant becau se the latter may influe nce the migration propertie s 
of the planet ( iPierens & Hurd 120051: IBaruteau & Massetll2008"l) . 
For global self-gravitating disks the treatmen t of the pot ential 
has been analyzed in more detail bv lHure & Pierensl(l2009l) . who 
calculated the required smoothing by comparing 3D and 2D disk 
models with specified density stratifications. For that purpose 
they compared the midplane value of the 3D potential to the 
2D case and estimated from this the required smoothing length. 
Additionally, they considered the whole extent of the disk for 
their integration. They found that e needs to be reduced for close 
separa tions ^ //, while for lo ng distances it approaches a finite 
value. iHure & PierensI (|2009|) give an extended list of smooth- 
ing prescriptions used in the literature and we refer the reader to 
their paper 

Here, we will reanalyze the required smoothing in 2D disk 
simulations. We show that the force to be used in 2D has to be 
obtained by performing suitable vertical averag es of the force . 
For the planet-disk case we extend the work by iMassetl (|2002|) 
and compare in detail the torque density to full 3D simula- 
tions. We present a method to approximately include the change 
in vertical stratification induced by the presence of the planet. 
With respect to self-gravitating disks we follow the work by 
iHure & PierensI (|2009|) and calculate the optimum smoothing 
length e by performing a vertical averaging procedure for the 
force between two disk elements. We show that this will be im- 
portant for fragmentation of gravitational unstable disks. 

In the next two sections we present the vertical averaging 
procedure and our unperturbed isothermal disk model. In Sect.|4] 
we analyze the potential of an embedded planet followed by the 
self-gravitating case. In Sect.|6]we summarize and conclude. 

2. The vertical averaging procedure 

Throughout the paper we assume that the disk lies in the z - 
plane and work in a cylindrical coordinate system {r,ip,z)- 
Starting from the full 3D hydrodynamic equations the vertically 
averaged equations, describing the disk evolution in the r - (p 
plane, are obtained by integrating over the vertical direction. The 
continuity equation then reads 



If-/ 



V ■ (pv) dz^Q, 



(1) 



where v and p are the 3D velocity and density, respectively. This 
is typically written as 



+ V. 
dt 



where 

2 = { pdz 

is the surface density of the disk and 



"4/' 



(pv) dz 



(3) 



(4) 



the vertically averaged 2D velocity in the disk's plane. The in- 
tegrals extend over all z ranging from -oo to +oo. The inte- 
grated momentum equation, considering only pressure and grav- 
itational forces, then reads 



2— = -VP 
dt 



-! 



pWdz, 



(5) 



with the vertically integrated pressure P - J p dz, where p is the 
3D pressure. From Eq. (|5) it is obvious that the change in the 
velocities is determined by the specific force (or acceleration) 
that is given by the ratio of force density (the integral in Eq.|5ll 
and surface density 2, i.e. 



/ = - 



/pV*Prfz 



(6) 



We here deal with the correct treatment of the last term in 
Eq. (|5]l, which describes the averaging of the potential that can 
be given, for example, by a central star, an embedded planet, or 
self-gravity 

^ = ^P. + ^Fp + ^Fsg . (7) 

From the derivation of the momentum equation (jSj it is clear 
that not the potential in the midplane matters, but a suitable av- 
eraging of the gravitational force. In the following we will an- 
alyze this averaging using density stratifications in the vertical 
isothermal case. Even though accretions disks are known not to 
be isothermal, we nevertheless prefer at this stage to use this 
assumption because it is still a commonly used approach that 
avoids solving a more complex energy equation. It leads to a 
Gaussian density profile that is also in more general cases a rea- 
sonable first-order approximation. Additionally, the irradiation 
of a central light source tends to make disks more isothermal in 
the vertical direction. More realistic structures will be treated in 
future work. 



3. The locally isothermal accretion disk model 

The vertical structure of accretion disks is determined through 
the hydrostatic balance in the direction perpendicular to the 
disk's midplane. In this section we present first an idealized 
situation, where the disk is non-self-gravitating and its struc- 
ture is not influenced by an embedded planet. This allows for 
analytic evaluation of integrals. Then we will consider more 
general cases. Using the thin-disk approximation and a gravi- 
tational potential generated by a point mass M* in the center, i.e. 
^, - -GMt/(r^ + z^y^^, the vertical hydrostatic equation for 
long distances r from the star can be written as 



]_dp_ 
p dz 



GM,z 

(r^+z^y- 



GMa 



^-nlz. 



(8) 



(2h) = 0, 



(2) 



To obtain the vertical stratification for the density, p, from this 
equation requires an equation of state p = p{p, T) and a detailed 
thermal balance, by considering, for example, internal (viscous) 
heating, stellar irradiation and radiative transport through the 
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disk. To simplify matters, often the so-called locally isothermal 
approximation is applied for the (numerical) study of embedded 
planets in disks or for self-gravitating disks. Using this approach, 
a costly solution of the energy equation is avoided by specifying 
a priori a fixed radial temperature distribution T{r). At each ra- 
dius, however, the temperature is assumed to be isothermal in 
the vertical direction. Using 



■■pel; 



(9) 



where the (isothermal) sound speed Cj is now independent of z, 
we obtain for the vertical disk structure a Gauss profile 



PG(r,(f,z) ^pQ(r,(f)e' 
where the vertical height H is defined as 



1 --2 

'2W 






(10) 



(11) 



and poir, (f) is the density in the midplane of the disk. For the 
vertically integrated surface density E one obtains in this case 



= I Poe -"'■ 



^ dz^ V2^//j 



Po. 



(12) 



In 2D, r - if, simulations of disks this surface density E(r, i^) is 
one of the evolved quantities. 

4. The potential of an embedded planet 

To illustrate the necessity of potential- (or force-) smoothing we 
consider in this section embedded planets in locally isothermal 
disks that are non-self-gravitating. First we will look at an unper- 
turbed disk whose structure is determined by the distance from 
the host star and is given by Eq. (fTOl i. Then we include the effects 
of an embedded planet on the disk structure. 

4.1. Unperturbed disk 

We consider the potential of a point like planet with mass Mp 
embedded in a 3D locally isothermal disk with a fixed Gaussian 
density profile, po- The planetary potential is given by 



*!'?('•) = -T 



GM„ 



GM„ 



{s^+z^Y 



(13) 



where r is the vector pointing to the location in the disk and r-g 
is the vector pointing to the planet. The vector s pointing from 
the disk element to the projected position of the planet, i.e. in the 
z = plane, is denoted by 



s = (rp -r)-(rp -r,e-)e-. 



(14) 



where e, is the unit vector in z direction and (■, ■) denotes the 
scalar product. 

The force acting on each disk element is then calculated from 
the gradient of the potential. Since it acts along the vector s, we 
can write for the modulus of the vertically averaged force density 



Fp(s) = - J p-^ dz = -GMpS J 



P 



{s^+z^)i 



(15) 



The force density (with units force per area) in Eq. (IT5i has to be 
evaluated at the centers of all individual grid cells, as illustrated 
in Fig. [U From there, either the total force of a disk element can 
be calculated by multiplying Fj,{s) with the disk element's area 




Fig. 1. Geometry of protoplanetary disk with embedded planet. 
We calculated the gravitational force from the planet (blue) on a 
vertical slice of the disk (red). For that purpose the gradient of 
the potential has to be vertically integrated along the dashed line 
that goes through the cell center To obtain the total force exerted 
on the shaded disk segment, this value has to be multipUed by its 
area in the r - tp plane. 



(rArA(p), or one computes the specific force /p = Fpfl., which 
can be used directly to update the velocities, see Eq. (|6]l. For the 
density we write 



pair, (fi, z) = po(r, (f) ■ Pz{z^/H^) , 



(16) 



where we have assumed that the vertical dependence of p is a 
function of the quantity z^/H^, as stated in Eq. ( fTOt . Substituting 
y - z/s in Eq. ( fTSl l and using a vertical stratification according 
to Eq. (fT6b . we find 



GMppo 



2Us) . 



(17) 



The dimensionless function Ip{s) is defined through an integral 
over half the disk 



Ip(s) = 

•I (1 + 



cY) 



(1+/)! 



dy, 



(18) 



where the normalized vertical distance y and the quantity c are 
given by 

and c = -^ . (19) 



y 



- and c = — . 
s H 

For the standard Gaussian vertical profile, i.e. Pzic^y^) 
exp(-icy ], the integral Ip(s) can be expressed as 



/n(^) 



1 2 

4^ ^''P\4 



^1 



^n 



(20) 



where K„(x) are the modified Bessel functions of the second 
kind. For illustrating purposes we present evaluations of the 
force correction function /p for simpler polynomial density strat- 
ifications in AppendixlBl 

In 2D numerical simulations of disks the above averaging 
procedure is typically not performed. Instead, an equivalent 2D 
potential is used in the momentum equation such that 



dt P ■ 



(21) 



We point out that a 2D potential with the property of equation 
Eq. (l2n i cannot be the result of an averaging procedure in general 
as in Eq. ^ because for realistic densities 



W: 



2D 



JpV'Vdz 



(22) 
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Fig. 2. Force correction function Ip{s) (solid line) resulting from 
an integration over the vertical structure of the disk, see Eqs. ( fTTl 
l20t . Additional curves indicate the corresponding function for 
the 2D potential (equivalent to Eq. (l24l i) using different values 
of the smoothing parameter fp. 

For ^Fp^ a simple smoothed version is often used, which reads 



m2D _ 

p 



GMr, 



(.2 + .2) 



(23) 



2V- 



where ep is the smoothing length to the otherwise point mass 
potential, introduced to avoid numerical problems at the location 
of the planet. We will refer to this functional form of ^i* as the e- 
potential, although it is sometimes named Plummer-potential as 
well. The force acting on each disk element is then calculated 
from the gradient of the potential. 

A finite fp regularises the potential and guarantees stable nu- 
merical evolution. Additionally, it serves to account for the ver- 
tical stratification of the disk. Comparing torques acting on a 
planet in 2D and 3D simulations it has been suggested for the 
Lindblad torques that e^ should be on the o rder o f the vertical 
disk height, specifically 0.7//, see iMassell (l2002h . He pointed 
out, however, that for the corotation torques a lower value may 
be appropriate. Hence, often a value of e = 0.6// is chosen. Here, 
we calculated the correct smoothing by a vertical average assum- 
ing an isothermal, vertically stratified disk. This will lead to a 
distance-dependent smoothing. 

In Fig.|2]we compare the force correction obtained from the 
vertically averaging procedure and the 2D smoothed e-potential. 
Specifically, we plot the function Ip{s) for the Gaussian density 
profile together with the corresponding function for the 2D po- 
tential, which reads 



I.{s) 



V2^// 



(^2 + 62)5 



(24) 



Because we are mainly interested in distances s up to a few //, 
we assume, that the disk height H does not change with radius. 
The unsmoothed ep = potential diverges as l/s for short dis- 
tances from the planet, leading to a l/s~ force as expected for 
a point mass. Since the value of Ip{s) remains finite for i — > 0, 
we see that a vertically extended disk reduces the divergence 
of the force to l/s. This is a clear indication that in 2D simula- 
tions the potential has to be smoothed for physical reasons alone, 
and that the assumption of a point mass potential will greatly 



overestimate the forces. As expected, the e-potential strongly re- 
duces the force for i — » and always yields regular conditions at 
the location of the planet. Additionally, it appears that the value 
ep - 0.6// overestimates the potential slightly for s > H (green 
curve in Fig.|2]i. 
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Fig. 3. Specific gravitational force Fp /E as a function of distance 
from the planet. Different approximations are shown: a) an ideal 
point mass (solid-red) that falls off as l/n; b) the exact averaged 
force according to Eqs. (l25l l and (l20t : c) the force according to 
the smoothed potential of Eq. ( l23l l using here ep = 0.7//; and d) 
a numerically integrated averaged force using 10 grid points in 
the vertical direction and a maximum Zmax - 5//. The force is 
normalized such that GM„ - 1 . 



Using the Gaussian isothermal density distribution, pc, and 
the surface density S from Eq. (fT2l i. the force density of the 
planet acting on the disk can be written as 



_GMpS2 



- Ipis) ■ 



(25) 



This expression could in principle be used directly in numeri- 
cal simulations of planet-disk interaction. However, even though 
/p(,s') can be solved in terms of Bessel-functions, in computa- 
tional hydrodynamics this evaluation is not very efficient, be- 
cause it has to be calculated once per timestep at each grid point. 
A possibility is to solve the integral Ip{s) in Eq. (fTSl l directly nu- 
merically using a limited number of vertical grid cells. In this 
case the integral is converted into a sum, where the exponential 
can be pre-calculated and tabulated at the corresponding nodes. 
For our purposes we found that only 10 vertical grid points give 
an adequate solution (see below). This is shown in Fig. [3] where 
we plot the specific force (acceleration) exerted by a planet on a 
disk element that is a distance s away. The force for a point mass 
falls off as l/s^ while for small radii the exact vertically averaged 
force shows a l/s behavior. Two approximations are displayed 
as well: The curve for the e- potential given in Eq. (l23l l. where 
we used a constant ep = 0.7//. The numerically averaged curve 
refers to a numerical integration of Ip{s) using 10 grid points in 
the z-range [0,5//]. Note that because of the finite discretiza- 
tion no additional smoothing is required. Increasing the number 
of grid points increases the agreement with the exact averaged 
force even more for shorter distances s. The scaling of the dis- 
tance with 1 /// and the normalization of the force make the plot 
independent of the used vertical thickness of the disk. While the 
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Fig. 4. Optimum value of ep in the smoothing potential of 
Eq. (l23l l as a function distance s from the planet. The shaded 
area illustrates which values of ep result in a force error of less 
than 10 % with respect to the exact averaged value. 



e-potential agrees well for s > H with the exact averaged force, 
our numerical approximation agrees to much shorter distances. 
In AppendixiBJwe show that simplified density distributions can 
lead to an equally good agreement. However, the advantage of 
the described vertical numerical integration procedure lies in its 
speed and in the fact that it leads directly to a regularized force 
for very short s. Also, it is easily generalizable, as will be shown 
in the next section. 

Because we can calculate I^{s) for the Gaussian profile nu- 
merically to any required accuracy, we can estimate the optimum 
e^is) value for the smoothed potential in Eq. ( l23l l for each point 
to obtain agreement with the exact averaged force. In Fig. |4]we 
display the correct e-p{s) value against the distance. The range of 
fp in which the smoothed potential produces an error of less than 
10 % in the force compared to the correct value is illustrated by 
the shaded region. Obviously for very short distances it is impor- 
tant to use the correct fp and for long distances the exact value 
of fp does not play an important role. Through a Taylor expan- 
sion for ,? — > oo of the denominator in Eq. (fTST l of I^{s) it can be 
shown that 



Um ep(s) - H , 



(26) 



which can be expected because the disk has a vertical extent on 
the order H. 



4.2. Taking the planet into account 

In the previous section we have analyzed the forces acting on the 
planet considering an unperturbed disk with a given scale height 
as determined by the star, see Eq. (flU . However, an embedded 
planet changes the disk structure, which will lead, for example, 
to a reduced thickness in the vicinity of the planet. This plays 
an important role for the torques acting on the planet. To esti- 
mate this effect we start from the vertical hydrostatic equation 
(O taking an additional planet into account. 



]_dp 
p dz 



GMa 



GM„ 



(r2+z2)5 (i2+^2)3 



(27) 



where s is again the distance from the planet. For a vertically 
isothermal disk this can be integrated 



Pp =poexp^---5— 



1 GM, 



■z' + 



GMp 

C'S 



(^2+^2)1/2 



1 



(28) 



where we assumed for the stellar contribution z «: r, as before. 
Using the previous vertical thickness H (Eq. (fTTT i) and the mass 
ratio q = Mp/M», this can be written as 



Pp-Poexp^-i^+^^ 



(^2+^2)1/2 



1 



(29) 



Because s is on the same order as z in the neighborhood of the 
planet, this equation cannot be simplified further To still obtain 
an estimate of the expected effects, we approximate the vertical 
density stratification by 



1 z^ kl 

Pn=Poexp<-|-— + — 



2//2 Hp 
where we define the reduced scale height near the planet 

4.2//2 



H„ 



qr^ 



(30) 



(31) 



This approximation for the vertical behavior of pp leads to the 
correct limits for the integrated surface density 2 in the limits 
for very short and long s and is a reasonable approximation in 
between. From the definition of //p in Eq. ( |3TI ) we find that the 
distance Si where the two scale heights are equal, 
given by 

r 2 \hl 



vl/2 



i.e. H„ - H, is 



(32) 



where h is the relative scale height, h — H/r, of the disk. The 
location St denotes approximately a transitional distance from 
the planet. For s < St the standard approximation of a Gaussian 
density distribution with the scale height H is no longer valid, 
and the influence of the planet dominates. 

Fig. |5] shows the exact (Eq. l29l l and the approximate ( l30t 
vertical density profiles for different distances to the planet and 
the comparison with the unperturbed Gaussian profile. For very 
short distances to the planet the effective height of the disk is 
much lower than in the unperturbed case. Obviously, Pp is only 
an approximation to pp, but it captures the change in thickness 
of the disk, as induced by the planet, very well. We will show 
below that using pjj in the force calculations yields a very good 
approximation to the exact case. 

To test how this new density-scaling influences the previous 
isothermal results for the planet's gravity, we calculated a cor- 
rected vertically averaged force using the modified density pp 
according to 



Fp{s) = -2GMpS f ^^—^dz. 

J (s^+z^)'- 



(33) 



This integral has to be calculated numerically using an approxi- 
mate Zniax- Depending on the model to be calculated (either 'ex- 
act' or approximate force), we substituted either pp from Eq. ( l29l l 
or the approximate Pp from Eq. ( l30l l. Since the presence of the 
planet alters the vertical height of the disk, Zniax depends on the 
distance s from the planet. To estimate Zm^x we first define an 
effective vertical scale height 



H, 



eff 



1 



1 



H'-^ H; 



-1/2 



(34) 



p/ 
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Fig. 5. Exact (see Eq. l29l l and approximate (see Eq. l30l l ver- 
tical density profiles for a disk with aspect ratio h = 0.05 for 
short distance i to a planet with mass ratio q - by. 10"^. The 
exact profiles are shown by the solid lines and the correspond- 
ing approximate ones in the same color but with dashed lines. 
For comparison the unperturbed Gaussian profile (see Eq. [10) is 
shown in black for this Hjr. 
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Fig. 7. Specific radial torque density in units of {dTldm\) (see 
Eq. [36b for embedded planet models using the e potential with 
various e^. The simulations use a given H - 0.05, and the .\:-axis 
refers to the radial coordinate where Op is the semi-major axis of 
the planet. The planet to star mass-ratio is (7 = 6 x 10"^. Other 
parameters of the simulations are stated in Sect. 14.31 
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Fig. 6. Specific force Fp/I. for different approximations. The 
first three curves are identical to Fig. [3]and are shown for com- 
parison. The 'exact' potential using pp is shown in yellow and 
the approximates solution, using the simplified expression for 
the density p^ and only 10 vertical grid points, is shown in pink. 
A tapering function near the center was used for the latter to 
regularize the force. 
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Fig. 8. Specific radial torque density in units of {dr/dm)o for 
embedded planet models for 2D models using e-potential (red 
curve) and a numerically integrated force F^ (see Eq. [25] l that 
takes into account a Gaussian vertical density distribution (blue 
curve). The green curve refers to a full 3D model using the same 
physical mo del parameter as t he 2D models. The 3D result is 
adapted from [Kiev et al.[ ( [20091) where the setup and numerics is 
described in more detail. 



which is an interpolation of the value at short and long distances 
from the planet. We choose to take Zm-^x = 6//eff in the 'exact' 
numerical evaluation where we use the density pp and 1000 grid 
cells. In contrast, we apply Zmax - 3//eff in the approximate nu- 
merical evaluation using p^ and only 10 vertical grid cells, see 
also Appendix [a] 

In Fig.[6|we display the results for the vertically integrated 
specific force F^/'L in various approximations. The first three 
curves (red, blue, green) correspond to those shown in Fig. [3]as 
an illustration. The 'exact' potential using pp is shown in yellow 
and the approximate solution using pjj in pink. The presence of 
the planet reduces the scale height of the disk, which leads to an 
enhancement of the force above the Gaussian, i.e., it diverges as 



s for short distances. To regularize the approximate force, we 
added a tapering function that reduces it to zero in the vicinity of 
the planet (pink dashed-dotted curve) such that it can be used in 
hydrodynamic simulations, see Appendix \A\ Clearly the range 
of applicability for short s is much improved over the simple e- 
potential. In contrast to the Gaussian approximation that uses the 
fixed thickness H, the force correction (with respect to the pure 
point mass potential) depends now on the mass of the planet as 
well, which enters through Hp. 
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4.3. Numerical simulations of planet-disk interactions 

To test the formulation of the force, we performed numerical 
simulations of embedded planets in two and three dimensions 
where we use an isothermal equation of state. For this purpose 
we solved the 2D and 3D hydrodynamic equati ons for a vis- 
cous g as. We used a se t up ver y similar to that of iKlev & Cridal 
(l2008l) and iKlev et alJ ( l2009l) . The planet with a mass ratios 
Mp/Mt - 6 X 10"^ is embedded at r = 1 in the disk with 
radial extent of 0.4 - 2.5. The disk is locally isothermal such 
that the aspect ratio H/r is constant. This implies a radial tem- 
perature gradient of T{r) oc r"', while in the vertical direction 
T is constant. As a consequence, the unperturbed vertical den- 
sity structure (in the 3D simulations) is Gaussian along the z- 
axis. The surface density falls off with radius as E oc r"'/^. We 
used a constant kinematic viscosity coefficient of v = 10"^ in 
dimensionless units. This setup is such that without the planet 
the disk is exactly in equilibrium and would not evolve with 
time. At the radial boundaries dampi ng boundary conditions 
were applied dde Val-Borro et al]l2006h . In our simulations we 
used H/r = 0.05 and evolved the disk after the planet's inser- 
tion for about 100 orbits. To test our improved treatment of the 
force, we varied the planet mass and the scale height of the disk 
in comparison models. For the standard parameter, h = 0.05 and 
q - 6x 10"^, we found for the transition distance Si ^ 0.35//. 

The embedded planet disturbs the disk and torques are ex- 
erted on it by the disk through gravitational back-reaction. These 
might lead to a change in the planet's orbital parameter The 
strength of these torques will depend on the applied smooth- 
ing of the gravitational force. To illustrate the effect, it is con- 
venient to study the radial torque distributi on per unit disk mass 
dr(r)/ dm, which we define here, following lD'Angelo & Lubowl 

mm . 



such that the total torque Ftot is given as 



27r 



f — 

J dm 



(r) 'L(r) rdr . 



(35) 



In other words, liF(r) is the torque exerted by a disk annulus 
of width dr located at the radius r and having the mass dm. As 
dY{f)ldm scales with the mass ratio squared and as {H/ry^, we 
rescale our results accordingly in units of 

where the index p denotes that the quantities are evaluated at the 
planet's position, with the semi-major axis Op. 

In Fig.Qwe plot dY{r)ldm obtained from 2D simulations us- 
ing the e-potential for the gravity of the planet and the standard 
fixed scale height H for the disk. The torque in this and the simi- 
lar following plots are scaled to {dr{r)/dm)o as given in Eq. ( l36l l. 
Results for five values of the smoothing length are presented. 
Obviously, the value of Cp has great impact on the amplitude of 
the torque density, and making the correct choice is important. 
We point out that the differences in the torque density also in- 
fluence the total torques that determine the important migration 
rate. Because Ftot consists of positive and negative contributions 
of similar magnitude, even small errors in dT{r)ldm can lead to 
large errors in F^t. For the range of e displayed in Fig.|7]we find 
a variation of Ftot larger than about a factor of 4. The best agree- 
ment in this case of H/r - 0.05 with 3D results is obtained for e 
in between 0.6 - 0.7//. 

In Fig. [8] we compare two different force treatments for 2D 
simulations to a full 3D run. Results for the e^ - 0.7// poten- 
tial are given by the red curve, and the blue curve corresponds 




Fig. 9. Specific radial torque density in units of {dT/dm)o for 
2D and 3D embedded planet models using different mass ratio 
and disk temperature. The first three curves represent models 
with the same disk temperature H/r - 0.05. The red curve cor- 
responds to our standard model, the blue curve has a reduced 
planet mass ^ = 3 x 10"^, and the third (green curve) corresponds 
to a full 3D run. The next two curves (yellow and purple) refer 
to models with H/r - 0.037. The 2D runs used our approximate 
density distribution p^ in eval uating the gravita tional force. The 
3D results were adapted from lKlev et al.l (l2009h where the setup 
and numerics is described in more detail. 



to the vertically averaged force F^ according to Eq. ( |25] |. which 
assumes a Gaussian vertical density profile. In magnitude the 
fp - 0.7// potential represents the vertically averaged results 
reasonably well, but it behaves differently close to the planet. 
The additional green curve correspon ds to a full 3D (l ocally) 
isothermal simulation as presented in iKlev et al] (|2009|) (their 
Fig. 10, purple curve). The 3D runs use an identical physical 
setup and a more realistic cubic -potential with a smoothing of 
Tsm = 0.5. Both 2D results are on the same order as the 3D run, 
but show small deviations that can lead to larger variations in 
the total torque, end hence the migration rate, because the posi- 
tive and negative contributions to the radial torque density are of 
similar magnitude. 

To test the applicability of our new procedure for treating 
the gravity in 2D simulations, we performed runs with different 
mass ratio and temperature in the disk. The results of 2D and 
3D simulations are shown in Fig. |9] where the torque density 
dr{r)/dm is plotted. The first set of simulations refers to H/r = 
0.05, where we compared the standard model in 2D and 3D to a 
run with half the planet mass. The next two curves show results 
of 2D and 3D simulations for H/r - 0.037. Despite the indicated 
differences in the parameter, all models used the same physical 
setup as described above. The 2D runs used our approximate 
density distribution pt, in ev aluating the gravita tional force. The 
3D results are adapted from iKlev et aD (l2009h where the setup 
and numerics is described in more detail. Firstly, all five curves 
show very similar overall behavior, confirming the scaling with 
{dr{r)/dm)(). The reduced amplitude of the H/r = 0.037 models 
is due to the onset of gap formation. Secondly, the agreement 
of the 2D and 3D runs is very good indeed. For example, upon 
varying the scale height, the change in shape of the curves is 
identical in 2D and 3D runs (yellow and purple) curve. We point 
out that the value of e to obtain the best agreement of the total 
torque Ftot in 2D and 3D simulations may depend on the value of 
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Fig. 10. Specific radial torque density in units of {dr/dm)o for 
2D and 3D embedded planet models using different mass ra- 
tio and disk temperature. The first two curves are identical to 
th ose in the previous figure, and the third is a fit presented 
bv iD' Angelo & Lubowl (l2010l) corresponding to the model with 
T{r) oc r"' and 2(r) oc r^'^-^. 




Fig. 11. Geometry of a protoplanetary disk for calculations with 
self-gravity. We calculate the gravitational force exerted by a 
vertical slice of the disk (blue) on another vertical slice of the 
disk (red), separated by a distance s. Two vertical integrations 
have to be performed along the dashed lines that go through the 
cell centers. To obtain the total force between the two segments, 
this value has to be multiplied by the corresponding areas, see 
also Eq. (3% . 



5. The potential of a self-gravitating disk 

Now we turn to full self-gravitating disks where a smoothing of 
the gravitational potential is required as well to account for the 
finite thickness of the disk. The potential at a point r generated 
by the whole self-gravitating disk is given by 



^sg(r) 



J \r-r'\ 



(37) 



Disk 



The smoothing required to obtain the p otential in the m i dplan e 
has been analyzed for this situation bv iHure & PierensI (l2009h . 
However, if one is interested in problems of fragmentation, it is 
more convenient to analyze the force between to individual ele- 
ments (segments) of the disk. Let us consider the force between 
two such disk segments that are a separated by the distance s, see 
Fig. [TT] The potential at the location r generated by a disk ele- 
ment located at r' which is a projected distance s away is given 
by 



^sg(r) 



/// 



Gpir',ip',z') 



dz'dA', 



(38) 



where dA' is the surface element of the disk in the r - (/?-plane at 
the point r' . The vector s is defined in analogy to Eq. (fT4l i and 
illustrated in Fig.fTTI 

The force at the location r generated by this vertically ex- 
tended disk element is calculated from the gradient of the poten- 
tial. The vertically averaged force density can then be written in 
analogy to Eq. ( fTST i as 



F..{s) 



/ 



P(r,(p,z)- 



^^Ps 



sg 



ds 



dz 



-Gs 



M 



p{r',(p',z')p(r,(p,z) 



{s'- + (z-z')^y- 



dz'dA'dz. (39) 



The evaluation of this integral depends on the vertical strat- 
ification of the density at both locations r and r' . As before, 
we consider locally isothermal disks. For weakly self-gravitating 
disks the density structure is then given by the Gaussian form in 
Eq. ([Tol l. However, similar to an embedded planet the disk's self- 
gravity will modify the vertical profile. Following our previous 
treatment of the embedded planet, we first analyzed the smooth- 
ing required for a disk that has an unperturbed Gauss profile and 
will then allow for modifications. 



H, because of the influence of the planet. Hence, it is always ad- 
visable to perform the simulations using the vertical integration 
of the force. 

To additionally validate our simulations we compare in 
Fig.[TO]our 2D an d 3D results obtained for the standard model to 
an analytic fit by ID' Angelo & Lubowl (1201 Ol) for the same disk 
parameter Although the fit has been developed for a smaller 
planet mass, the agreement of the two 3D results is excellent. 
This is interesting because our planet mass of 20 Meaith is al- 
ready in the range where non-linear effects should set in. The 2D 
torque shows the same amplitude, but small differences are visi- 
ble just inside the planet. In isothermal disks the flow close to the 
planet is not in exact hydrostatic equilibrium any more. However, 
full 3 D high-resolution isothermal simulations (ID' Angelo et alJ 
l2003h have shown that vertical velocities are only large inside the 
Roche region of the planet. This explains the good agreement of 
the 2D approximation with the full 3D case. 



5.1. Unperturbed disk 

For a locally isothermal disk and with Eq. ([16} we obtain 

Psg(s) = -Gpoir, ip) I I poir', ip') dA' ■ 2h^{s) , 

where we defined the function hg{s) by 



(40) 



2 J J (I + (> 



y )p-ic y ) 



— oo —oo 



il + {y-y')^r- 



— dy' dy , 



(41) 



where the normalized vertical distance y and the quantity c are 
given by Eq. ( fT9b again assuming a constant H. This integral 
cannot be calculated directly for the standard Gaussian profile, 
and we will evaluate it numerically. 
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Fig. 12. Force correction function Isg{s) multiplied by s result- 
ing from an integration over the vertical structure of the di sk, see 
Eqs. ( 140141b . Additional curves indicate the corresponding func- 
tion for the 2D potential (Eq. l24l i using different values of the 
smoothing parameter fsg- 

In 2D numerical simulations usually a simple smoothed po- 
tential is used instead of calculating the correct averaging. This 



-potential reads as 



n"(*) 



-// 



GS(rO 



dA' , 



(42) 



where ejg is the smoothing length. The force acting on each disk 
element is then calculated from the gradient of ^l^. 

In Fig.[T2]we compare the force correction obtained from the 
vertically averaging procedure and the 2D smoothed potential 
for a disk height that does not change with radius. Because hgis) 
diverges for short distances, we multiply it by s. Because we are 
interested in local effects with distances i up to a few H, we as- 
sumed for the plot that the disk height H is constant. Comparing 
this result to the corresponding force correction function for the 
planet in Fig. |2l it is obvious that for the self-gravitating case a 
larger smoothing is required for the e-potential than in the planet 
case. Here, a value of e = 1.2 seems to be a good choice. Lower 
values, already e = 0.6, considerably overestimate the force. We 
attribute this larger required smoothing with the double vertical 
averaging that has to be performed in this case. 

From the numerically calculated hgis) for the Gaussian pro- 
file we can calculate the best esg(i) value for the smoothing po- 
tential in Eq. ( l42l i at each distance s to obtain the right force cor- 
rection value. In Fig.[T3]we display the optimum e^gis) value ver- 
sus distance. The range of esg over which the smoothing potential 
produces an error of less than 10 % is shaded. Cleary for short 
distances it is crucial to use the correct value of esg, whereas for 
long distances the influence of esg becomes negligible. Because 
we now have to account twice for the vertical extent of the disk 
here compared to the planet case, we obtain a higher limiting 
value of 

limesgC,?) = V2//. (43) 

This value was obtained through numerical calculation up to 
the sixth significant digit for s = 1000//, and we verified it 
through a Taylor expansion of the denominator in Eq. (l42t . As 
before, it is interesting that the required optimum smoothing re- 
mains finite, even for v ery long distances. This result agrees with 
lHure&Pierensl(l2009h . 




Fig. 13. Optimum value of esg in the smoothing potential of 
Eq. (l42l i as a function of distance s using a Gaussian vertical 
stratification. The colored area illustrates which values of esg re- 
sult in an error of the force correction value of less than 10 % . 



5.2. Taking the disk into account 

Now we consider the correction for a disk where self-gravity 
has modified the vertical structure. For that purpose, we con- 
sidered the case of a pure self-gravitating disk, neglecting the 
gravitational potential of the central mass. Assuming a large disk 
with a slowly varying surface density £(r), then the derivate with 
respect to z of the disk potential given in Eq. ( |37] | simplifies 
(,Mestel,1963,) to 

—^ = 27rGS(r) , (44) 

oz 

and consequently the hydrostatic equation dH) changes to 



]_dp 
p dz 



-InGYir) . 



(45) 



For a vertically isothermal disk this can be integrated (ISpitzeJ 
[T94l to 

(46) 



(47) 



Psg(r, (fi, z) = pair, ip) cosh I ^ h 
where the vertical scale height H^g is defined by 



//sir — 



c? 



Csi^K 



ttGI, nGI. 



H=QH, 



with the Toomre parameter Q (lToomrelll964l) . Fig. [T4l shows 
the changed vertical density profiles in the self-gravitating case 
compared to the unperturbed Gaussian for H = H^g or Q = 1, 
respectively. The self-gravitating profile is steeper and there- 
fore, for equal surface density, the mass is consequently located 
nearer to the midplane of the disk. This is expected because 
the vertical component of the disk's gravitational potential is, 
with InCLir) ~ GM/r^, by a factor of about r/H larger than in 
Eq. (|8]i, and so the mass is more concentrated toward the mid- 
plane. 

Now we can calculate the force correction function hg{s) 
of equation (HTl i with the self-gravitating vertical density pro- 
file psg. In analogy to the unperturbed case, this can be used 
to calculate an optimum esg for the e-potential. In Fig. [15] we 
display the correct e^gis) value against distance. The optimum 
esg for this value of Q is always lower than our previous un- 
perturbed Gaussian case. For the limit s — » oo we find a value 
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Fig. 15. Optimum value of fsg in the smoothing potential of 
Eq. (l42l i versus distance s using a vertical stratification caused 
by self-gravity with Q - \. The colored area illustrates which 
values of fsg result in an error of the force correction value of 
less than 10%. 



about 10% lower. This is in consistency with the lower effec- 
tive vertical scale height, because more mass is located near the 
midplane. 

In most self-gravitating disks the vertical structure will be 
affect by both the central mass object and the self-gravity of the 
disk. Then it is to be expected that the correct fsg is a value be- 
tween the two extreme cases. The combined situation where self- 
gravity and the central mass both contribute has been considered 
by ILodatOi (i2007i) . For clarity, we treat the two cases separately 
here. 

5.3. Numerical simulations of self-gravitating disks 

To demonstrate the influence of tsg in numerical simulations, we 
analyzed the effects of different values of esg on the fragmenta- 
tion conclusions of a self-gravitating disk. For that purpose, we 
adopted a disk model fromL Bamteau et al. ( 201 1) . Table[T]shows 
all important disk parameters. We sim ulated the disk twice using 
the ADSG version of the FARGO code dBaruteau & MassetlllOOa 
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Fig. 14. Vertical density profiles of the disk with for a self- ^ 
gravitating disk (see Eq. |46] | with Q - \. For comparison the g lo' 
unperturbed Gaussian profile (see Eq. [TOt is shown in red. 
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Fig. 16. Radial dependency of surface density {top), midplane 
temperature (middle) and Toomre parameter (bottom) of the two 
disk models with fjg = 0.6// and esg - 0.006// for different 
timestamps. At r = a both simulations match because they have 
identical start conditions. In the first 500 years they behave very 
similar but then their further evolution diverges drastically. 



lM.asset"2000^ for 5000 years with values of 0.6// and 0.006// for 
esg. Both models include the self-gravity of the disk and a simple 



Table 1. Parameters of the disk model. The top entry refers to 
host star The following give the disk properties, and then the 
initial disk setup and the computational parameter are given. 



Star mass (Mp,.,n,ary) 


1.0 Mq 


Disk mass (M^isk) 


0.4 Mq 


Adiabatic index (y) 


5/3 


Mean-molecular weight Qi) 


2.4 


/3-Cooling (/?) 


20 


Initial density profile (S) 


oc r 


Initial temperature profile (T) 


oc r~ 


Initial disk aspect ratio (H/r) 


0.1 


Grid (A^, X A'^) 


512 X 1536 


Computational domain (/?,„!„ -/?max) 


20-250AU 



10 
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jS-cooling model, which is defined by 



dt 



EQ. 



(48) 



where E is the internal energy, Q the angular velocity and /? a 
constant. The disk must cool fast enough to be able to frag- 
ment, wh i ch is otherwise prevented by compressional heating. 
Gammi d (|200lh showed that this is the case for /? < 3 and 
Rice et d] (l2005h found later a dependency from the equation 
of state for and suggest ed a yS < 7 for y = 5/3. 

iMeru & Batd (1201 lah pointed out that the previous results 
had not converged with increasing resolution and that the crit- 
ical value, yScrit, may be higher than previously thought. They 
measured an /?crit of ~ 18 for their highest resolution. Because 
we use yS = 20 in our models, we do not expect them to fragment 
in any case. 

Another stability criterion can be described by the Toomre 
stability parameter Q, which is defined by 



Q 



ttGI. 



(49) 



where C s is the sound s peed in the disk and k the epicyclic fre- 
quency. iToomrd ( 1 19641) showed that for values of 2 ^ 1 an ax- 
isymmetric disk should be stable. As shown in Fig. [16] the mod- 
eled disks have an initial Q value of 1.85 - 2.2 depending on the 
radial distance to the star, and thus should not be fragmenting. 

In the first 500 years both disks cool down in the central parts 
of the disk from about 120 K to about 30 K, resulting in Q values 
of 1 at the inner edge of the disk, which means that they are not 
Toomre-stable anymore, but still should not fragment because 
the cooling constant /3 is higher than the critical value required 
for fragmentation. After 500 years both simulations start to dif- 
fer The fsg = 0.006// model fragments within about 500 years 
in the inner region of the disk whereas, the e^g = 0.6// model 
needs about 1000 years to start developing small spiral arms in 
the inner region. After 5000 years two fragments have survived 
in the e^g = 0.006// disk with fragment masses of 16.5 Mj^p and 
8.3 Mjup. The e^g = 0.6// model only shows spiral arms and no 
signs of fragmentation. Even after 50000 years we did not ob- 
serve any fragments. Fig. [T7]shows the surface density distribu- 
tion of both models after 5000 years. 

In Section 15.11 we suggested an e^g on the order of unity to 
obtain results that can compare to 3D simulations. Because the 
esg = 0.6// model did not fragment as predicted by the stability 
criteria, this seems to support the validity of our estimate for e^g. 
The very short esg in the tsg - 0.006// model overestimates the 
gravitational forces on short distances (see Fig. [12) and therefore 
excites disk fragmentation. 



6. Summary and conclusions 

We analyzed the smoothing of gravity in thin 2D disk simula- 
tions for the embedded planet and self-gravitating case. Starting 
from the vertically averaged hydrodynamic equations, we first 
showed that the gravitational force has to be calculated using a 
density-weighted average of the 3D force, see Eq. ([5). Because 
this depends on the density distribution, there cannot be a gen- 
eral equivalent 2D version of the potential (Eq. l22l i. To be able 
to explicitly calculate the averaged force, we first used a locally 
isothermal disk structure in which the vertical density stratifica- 
tion is Gaussian. 

For the embedded planet case the resulting force can be cal- 
culated analytically. In full 2D hydrodynamic simulations we 



compared the resulting torque density acting on the planet for 
the e-potential ( 123] ) and the 'exact' averaged force. We found 
that the overall magnitude of the torque is best modeled using 
a smoothing of e = 0.7//, while there remain significant dif- 
ferences to the full 3D case, also in the total torque. Taking the 
modification of the density stratification induced by the planet 
into account leads to a much reduced vertical thickness of the 
disk in the vicinity of the planet. We presented a simplified an- 
alytical form for the modified vertical density stratification with 
a planet, the details of the numerical implementation are given 
in the appendix. Using this in 2D simulations leads to very good 
agreement of the torque density with full 3D calculations of em- 
bedded planets on circular orbits. By varying the disk height and 
the mass of the planet, we showed that the torque density scales 
as expected with {dr{r)/dm)o, see Eq. ( l36b . 

Because the modified density approximation that includes 
the planet is based on vertical hydrostatic equilibrium, it is not 
clear however, whether it is valid for planets on non-circular 
or inclined orbits. Here, the variations occur on the orbital 
timescale and to establish hydrostatic equilibrium, the thermal 
timescale must be on the same order In the situation of multi- 
ple planets that may interact strongly, the same restrictions may 
apply. Despite these restrictions, we believe that using this pre- 
scription will enhance the accuracy of 2D simulations consid- 
erably. We expect that our procedure can be generalized to the 
radiative case using suitable vertical averages, but this needs to 
be developed. 

For the self-gravitating case we showed that the required 
smoothing, e a: //, is even larger than in the planetary case. 
We attribute this to the necessity of a double averaging over 
the vertical height of the disk. Owing to the complex integra- 
tion, the integrals cannot be solved analytically in this case. 
Taking into account self-gravity lowers the required smoothing 
because the vertical scale height is reduced due to the additional 
gravity. In more strongly self-gravitating systems, which have a 
Toomre parameter Q ~ I, non-axisymmetric features may occur 
Because the standard self-gravity solvers require a smoothing 
that scales with radius, one has to take an approximate average 
in non-axisymmetric situations. The same applies for disks that 
are close to the fragmentation limit. As shown by our last ex- 
ample, the choice of smoothing may affect the conclusions on 
whether the disk will fragment or not. Through detailed com- 
parisons with full 3D simulations a suitable smoothing may be 
found. 

Appendix A: Integration of the force density 

Here, we briefly outline a numerically fast and convenient 
method to vertically integrate the force for the embedded planet 
case. Specifically, we plan to evaluate the force density Fj, in 
Eq. ( l33b , which reads 






(A.1) 



As pointed out in the text, we used for the (approximate) numer- 
ical integration a maximum z of Zmax - 3//eff , with H^ff given by 
Eq. ( l34b . The interval [0,Zniax] is divided into A^- equal intervals 
with the size Az = Zmax/N-. The integral in Eq. dA.lb is replaced 
with the following sum 



/ 



Ppdz 

{s'-+z^)^- 






Ppjzk) Az 



(A.2) 
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Fig. 17. Surface density of two disks after 5000 years which started from the same initial condition but with different values for the 
smoothing parameter esg of the gravitational potential. The gray circle shows the computational domain. 



where the A^^ nodes are located at z^ = (k - I /2)Az. We intro- 
duce a small smoothing, r,, here to keep the sum regular at short 
distances s. In the simulation presented in the paper we used 
Ts - 0.1/?Hiii throughout, which is much shorter than the unper- 
turbed vertical height H. In the 2D hydrodynamic simulations 
we used A^; = 10, which makes the method feasible numerically. 
The specific force, or acceleration (Fp/2), is then obtained 
by dividing through the integrated surface density 



Appendix B: Approximate vertical density profiles 



2 = 2^ppfo)Az, 



(A.3) 



k=i 



using the same nodes. In 2D simulations S is one of the evolved 
quantities and this relation can be used to calculate the other- 
wise unknown midplane density po- Because the vertical num- 
ber of grid points is very small (A^; - 10), the method is reason- 
ably fast and can be used in numerical planet-disk simulations. 
Additionally, it only needs to be evaluated in the vicinity of the 
planet and could be omitted farther out. Nevertheless, despite 
the coarseness of the integration, the agreement with the 'exact' 
force and torque is excellent. For numerical stability we smooth 
the resulting force using the following tapering function 



/taper I -^ J — 



1 



exp[-(,9-rt)/(0.2r,)] + I 



(A.4) 



with the tapering cutoff-length r^, for which we used in our sim- 
ulations rt - O.lRum- The purple curve in Fig. |6]exactly corre- 
sponds to the procedure described here, using the stated param- 
eter. 

Finally, we point out that for the stellar contribution to the 
potential 

„.„ = _ii^ = -^i^, (A.5) 

If- '■•I (r'-+z-r- 

a very similar vertical averaging needs to be performed. 
However, because the vertical profile is always Gaussian, the 
nodes can be precomputed and the exponential function has only 
to be evaluated once for the A^^ nodes per grid point and timestep. 
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Fig. B.l. Gaussian vertical profile p. (solid line) compared with 
the parabolic vertical profile p-' and fourth-order vertical profile 
p~ . They all have the same area below the curves and yield the 
identical surface density. 



To simplify some estimates and obtain an idea of the func- 
tional behavior of the forces, it is useful to study simpler vertical 
density stratifications. Here, we present results for a parabolic 
and quartic behavior The corresponding density stratifications 
read 



pf 



for the parabolic form, and 



1 



-f— f 



^(4) 



^ 2l//(4)J "*" 16\HW) 



(B.l) 



(B.2) 
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for the quartic form. The vertical heights //'^* and //*"*' of the 
models are specified such that the corresponding surface and 
midplane densities match the isothermal case, see Eq. (fT2l i. We 
obtain //<2) ^ 3/4 V^// and //<4) = 15/16 y/^H. Fig. lO 
shows all three vertical density profiles in comparison. 
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Fig. B.2. Force correction functions for different vertical den- 
sity profiles. /p(s) is the numerically calculated force correction 



(2), 



M), 



function for the Gaussian vertical profile p,. /J, (s) and /p (s) 
are the analytically calculated force correction functions for the 

(2) 

parabolic vertical profile p. and the fourth-order vertical profile 

(4) 
Pz ■ 
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These density stratifications are then used to calculate the 
vertically integrated force and to obtain the corresponding force 
correction functions. Here, the vertical integrations extend to 
that Zmax value where the density vanishes. For the two distribu- 
tions we find zJnax = V2//'^^ and Zmlx - '2H^'^\ For the second- 
order integral we find 



e\s) 



1 



1 



2 ' 2 
and for the fourth-order integral 



Ct arcsinh 



(I) 



P ^' 16 



V^ 



+ 4 



1 + 



16 



arcsinh - 



Q 



(B.3) 



(B.4) 



Here, we defined C2 - s/H^^^ and C4 = s/H^'^\ As Fig. lB.2l illus- 
trates, for the unperturbed disk (without a planet) these functions 
agree reasonably well with the Gaussian value also for smaller 
separations s. For long distances all studied force correction 



f(2) 



M) 



functions (/p,4,/p ' and /p ) approach each other These sim- 
pler profiles may be also useful in the study of self-gravitating 
disks. 
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